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ABSTRACT 



Aims. Determination of horizontal velocity fields on the solar surface is crucial for understanding the dynamics of structures 
like mesogranulation or supergranulation or simply the distribution of magnetic fields. 

Methods. We pursue here the development of a method called CST for coherent structure tracking, which determines the 
horizontal motion of granules in the field of view. 

Results. We first devise a generalization of Strous method for the segmentation of images and show that when segmentation 
follows the shape of granules more closely, granule tracking is less effective for large granules because of increased sensitivity 
to granule fragmentation. We then introduce the multi-resolution analysis on the velocity field, based on Daubechies wavelets, 
which provides a view of this field on different scales. An algorithm for computing the field derivatives, like the horizontal 
divergence and the vertical vorticity, is also devised. The effects from the lack of data or from terrestrial atmospheric distortion 
of the images are also briefly discussed. 
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1. Introduction 

Determining the flows on the surface of the Sun has trig- 
gered many efforts over the last decade. Of special con- 
cern was determination of granules motions that can re- 
veal horizontal flows on scales larger than, typically, twice 
their horizontal size ~ 2000 km. This scale is small enough 
to provide observers of the Sun's surface with a detailed 
sampling of the large-scale flows, such as supergranula- 
tion, and therefore makes it very interesting to determine 
granule motions. 

However, from the viewpoint of fluid mechanics, gran- 
ules are not passive scalars whose motions trace that of 
the fluid; rather, they are structures in intensity or, assum- 
ing perfect correlation with temperature, coherent struc- 
tures of temperature. But the evolution of the tempera- 
ture field is the result of radiative, as well as advective, 
processes. It is only in case that the latter dominates that 
the motions of granule can be associated with horizontal 
flows. Using numerical simulations in a large horizontal 



box, we have shown that granule motions are highly cor- 
related wit h horizontal flows wh en the scale is larger than 
^^2500 km ( Rieutord et al.|[200lh : below this scale, granule 
motions should be considered as (solar) turbulent noise. 

Once the equivalence of plasma motion and granule 
motion is assumed, one is left with the problem of mea- 
suring the latter motion. This is not an easy task ow- 
ing to the small angular size 1.3") of the structures. 
Ground-based observations are sensitive to atmospheric 
turbulence, while space observations are expensive owing 
to the (relatively) large aperture needed for resolving gran- 
ules. 

Basically, two techniques have been used to mea- 
sure horizontal velo city fie l ds: e ither the tracking of 
individua l granules (I Stroud 1 19951) or local correlation 
tracking ([November fc SimonI 19881. Th e results of these 



two technique s have been compared (Simon et al.l 11995 
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1999[ ) and found to bro adly agree; in 
the test using numerical simulations (jRieutord et al. 
20011) . they show the same degree of correlation with 
the actual plasma flows. However, detailed examina- 
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tions ( Simon et al.ll995t Roudier et al .11 19991 ) have demon- 
strated worrying differences, especially when field deriva- 
tives like vertical vorticity and divergence are computed. 

In fact from the point of view of signal processing, 
these two methods differ fundamentally. On the one hand, 
granule tracking emphasises the importance of the gran- 
ule, gives no signal in between granules, and yields a ve- 
locity field that is sampled randomly following the distri- 
bution of granules. On the other hand, local correlation 
tracking (LCT hereafter) treats granules and intergran- 
ules on an equal footing and yields a velocity field on a 
regular grid. Broadly speaking, the two methods differ in 
the interpolation process, which unfortunately influences 
the final result. 

In this paper we present and analyse in some detail an 
algorithm based on granule tracking which is able to give a 
reconstruction of the velocity field at all scales larger than 
the sampli ng scale. This algorit hm has already been in- 
troduced in! Roudier et al ] (|l999l ) in a preliminary version. 
We call it CST for coherent structure tracking to under- 
line its relation with the physics lying behind it. Such an 
algorithm is close in its principles to particle-imaging ve- 
loci metry (PIV), as used in experimental fluid mechanics 
(e.g. [Adria3 [2003). 

We developed this algorithm for three reasons: the first 
is obviously because it gives a different view of the data 
than does the LCT algorithm, since many interpolation 
probler ns may influence the final results (see the discus- 



sion m iPotts et al.l 120031 ) . The second one is that it may 
be used on raw data and gives an estimate of the error 
introd uced by atmospheric distortion (see the companion 
paper, iTkaczuk et al. 2007L in which this point is devel- 
oped). Finally, it offers the possibility of selecting specific 
structures according to their nature, size, lifetime, etc. and 
of studying their motion. 

In the next section we discuss the different steps of the 
algorithm, especially the segmentation and interpolation 
processes, and also point out the effects of regions lacking 
data. Discussion and conclusions follow. 



2. The CST algorithm 

Before describing the different steps in some detail, let us 
recall the five main steps of this algorithm: 

— segmentation of the image and granule identification 

— measurement of velocities at granule locations 

— reconstruction of the velocity field 

— calculation of field derivatives (like the z-component of 
the vorticity and divergence) 

— estimation of the noise. 

We now go into detail for each one in turn. 

2.1. Segmentation and granule identification 

To identify a granule one needs a criterion with which to 
decide whether a given pixel belongs to a given structure 
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Fig. 1. Comparison of the result of different segmenta- 
tions on the same test image. From top to bottom and 
left to right we have : the original image, the result of 
Strous's method, the result of the BW method, the re- 
sult of the proposed method. Note that BW leaves some 
granules undetected. 



or not. This criterion needs to be local in order to avoid 
threshold effects due to large-scale variations in intensity, 
which either come from terrestrial atmospheric effects, so- 
lar acoustic waves, or even magnetic fields. 

A classical criterion is based on detecting local max- 
ima of the inten sity through the curvature C = h+i — 
h - {h - h-i) (|Strouslll995l I Roudier et aTlll999l ). This 
criterion has the advantage of being simple, robust, and 
therefore quite efficient. However, comparing the detected 
patterns with the original image shows that this criterion 
underestimates the size of the granules. It is therefore in- 
teresting to know whether this criterion can be improved. 
An objective test of this improvement will be that the 
lifetime of the granules is increased. 

Another ni e thod has been proposed by 
Bovelet fc Wiehr ( 200l[ ) (hereafter referred to as BW) 
with a multiple-level tracking algorithm. It is based 
on the use of multiple-threshold levels applied to the 
intensity image. The granules detected at a high level 
are gradually extended to adjacent pixels whose inten- 
sity exceeds the lower level, while keeping a minimum 
distance with respect to other granules. This approach, 
which is very sim i lar to a watershed- based segmentation 



(jVincent fc Soilld Il99ll : ISoilld Il999[ ). yields sizes and 



shapes of granules that conform more to direct observa- 
tion of the image. Nevertheless, the number of detected 
granules is less compared to the Strous algorithm because 
this method is based on the intensity. 
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Fig. 2. Comparison of the effect of different segmenta- 
tions on lifetimes of granules. BW refers to the Bovelet 
and Wiehr method while CD refers to our method (see 
text); the associated number is the extension threshold 
parameter text- 



The Strous-curvature criterion is more efficient at sep- 
arating the granules. It consists in selecting the pixels 
whose local minimal curvature is not negative. This cur- 
vature is calculated in the four directions defined by the 8 
neighbours of the considered pixel. The underestimation 
of the granule size stems from real granule's extension be- 
ing greater than the positive curvature region. We thus 
propose a new segmentation algorithm (hereafter termed 
CD) that combines both ideas of BW and the Strous al- 
gorithm. It consists in the following steps: 

— Calculation of the "minimal curvature image" : for 
each pixel, the minimal curvature among the four di- 
rections is computed. 

— Detection of the granules as non negative curvature 
pixels in the minimal curvature image. 

— Extension of the detected granules with points whose 
minimal curvature value is greater than a given (nega- 
tive) threshold t^xt^ while keeping a minimal distance 
of one pixel between each pair of granules. 

The last step can be reached using the watershed algo- 
rithm on the minimal curvature image, with an additional 
condition requiring that the curvature remains above t^xt ■ 
This new approach leads to a segmentation with the 
same granules as Strous's approach, but with a control- 
lable size. Strous's segmentation is obtained with t^xt = 0. 
Decreasing the value of text extends the granules. In Fig.[T] 
we illustrate the discussed segmentation method using 
Pic-du-Midi dat£0. This figure illustrates the way our seg- 
mentation extends that of Strous and closely follows the 
shape of granules. 



To study the influence of the segmentation on the life- 
time of granules, we plotted the statistic of the lifetime 
for the three methods : the Strous method, BW method, 
and our proposed method. For our method, we took three 
different threshold values: text = 0, text = —0.1, and 
text = -0.3 (Fig. El). 

This figure shows that, broadly speaking, the segmen- 
tation does not influence the statistics of lifetimes very 
strongly. However, some differences can be noticed in the 
detail. The BW algorithm, by detecting less granules, 
shows a deficit of short-lived and long-lived granules; on 
the other hand, our algorithm eliminates long-lived struc- 
tures when used with too low a threshold. We understand 
this behaviour as the result of enhanced splitting of larger 
structures. 

In conclusion, Strous's algorithm seems the most effi- 
cient for our purpose, and it can be improved in the way 
we described, but at the price of increasing the computa- 
tion effort a lot. 

Once the image has been segmented, each granule 
needs to be identified. This operation, although very sim- 
ple, can be quite time-consuming since all pixels should be 
tested at least once. The most efficient way we have found 
to deal with this operation is to use a recursive algorithm, 
letting granules grow from a single pixel. A pixel belongs 
to a granule if it shares at least one side with another pixel 
of the granule. 

2.2. Measuring the velocities 

Once the granules have been identified, the (x, y) coordi- 
nates of their barycentre are computed. Hence, each im- 
age is converted into a set of points Xi^n describing the 
position of the granules at time These data may be 
completed by the set of granules surfaces, shapes, etc. 

The set of points {Xi_„} is then divided into trajecto- 
ries 



{^»(fc),"}„. 



<n<ni 



This notation means that granules of index i{k) are in fact 
the same granule as the one that follows the k^^ trajectory; 
it appears at time t{ni) and disappears at time t{ni). 

Trajectories are identified by comparing each posi- 
tion Xi on two consecutive images and putting nearest 
neighbours together provided their position does not dif- 
fer more than a given threshold that is determined by an 
upper bound on the velocity. Typically, we reject velocities 
higher than 5 km/s. 

If one disposes of a long time series of images, it may 
be useful to determine the time evolution of the velocity 
field. For this purpose a time window of width At needs to 
be used and trajectories are restricted to time windows. 
Hence, for a given time window, one derives a set of N 
trajectories: 



As for all examples needing solar images, we used the series 
obtained at Pic-du-Midi on 20 September 1988. 



<n<ni 



t< t{n,), t{nt) < t + At\ 
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The values of the velocities are of course not uni- 
formly distributed in the field of view and we need to 
know how they constrain the velocity field at a given res- 
olution: small-scale components are weakly constrained, 
while large-scale ones are highly constrained. The maxi- 
mum resolution for the velocity field is given by the den- 
sity of trajectories and can be estimated by the mean dis- 
tance < d > between the Xk- As the time window At 
is increased, the maximum resolution increases according 
to the law < d >cx l/\/At as clearly shown by Fig. |4l 
This law arises because granules cover the Sun's surface 
permanently. From it, we can derive the maximal spatial 
resolution for a given time resolution. Indeed, from Fig. HI 
it turns out that the mean volume of space-time occupied 
by one granule is ~ 1200 Mm^s. 



42.6 
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Fig. 3. Trajectories of granules in a five-minute time se- 
quence. Top: Granule displacements in a small subfield. 
Bottom: enlarged view of a granule trajectory. Data ar e 
from the 1988 Pic-du-Midi series fe.g. iMuller et al.lll992l ). 



existing in the field during the time [t, t + At]. We show 
an example of such a set of trajectories in Fig. [31 We 
clearly see from the enlarged view that granule motion 
is dominated by an erratic motion that mixes (Earth) at- 
mospheric noise with the turbulent random flow of solar 
convection. 

From this set we derive a mean velocity associated with 
each trajectory; the k^^ trajectory gives the velocity 



t{n2) ~ t{ni) 



which we associate with the mean position of the trajec- 
tory 



1 



"2 

^^^^^^ ^ 

n— ni 



Hence we end up with the set 

which describes the velocity field during the time interval 
[t,t + At]. 



10 

At (min) 



Fig. 4. Mean distance between velocity vectors as a func- 
tion of the time window used for the measurement. The 
line shows the law < d >cx l/y/At. 



This value is useful for determining the highest time 
resolution that can be allowed fo r the large-scale veloc - 
ity fields. Indeed, as pointed out in ' Rieutord et al.' (2001), 
granules cannot be used to trace plasma flow under a scale 
of 2.5 Mm (except in the case of very rapid flows like "ex- 
plosion" of granules); thus the determination of a large- 
scale flow needs a mesh size not smaller than 1.25 Mm. 
This means that each grid point will have a have a trajec- 
tory only if At > 768s. Conservatively, we estimate that 
the time resolution cannot be higher than 15 min. 

Often time resolution is not needed and therefore At is 
much larger than 15 min. In such a case, several velocities 
may be given by granules appearing at a given place on the 
grid. We then use the average velocity of the granules as 
the measure of the local velocity. This introduces another 
quantity, namely the rms fluctuation around this mean. 
This rms velocity is a measure of the proper velocity of 
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granules and therefore a proxy of the local strength of 
turbulent convection. 



2.3. Derivation of the velocity field by MRA 

Once the velocities in the field of view are determined, 
we need to know what kind of flow they represent: for in- 
stance, vortices, shear layers, diverging sources, etc. For 
this purpose we need to determine the best continuous 
differentiable fleld that approximates the data. The de- 
termination of such a field can be done in various ways, 
but we wish to introduce no information or a minimum of 
into this process. We also wish to avoid the propagation 
of errors and noise in the field of view. 

We have found that wavelet multi-resolution analysis 
(MRA) is an interesting tool for this purpose, because it 
gives a decomposition of the signal at all the scales allowed 



by the size of the box and therefore gives a g ood view of a ll 
the components of the signal (Mallat 198a: lMeveiill987 ). 



The basic idea of a "multi-resolution representation" 
of L^(IR) is to project a signal / on a "wavelet orthonor- 
mal basis" of L^(5R), at which point it is possible to ex- 
tract the "difference in information" between two succes- 
sive approximations of the signal (approximations at the 
resolution 2^ and 2-'"'"^). The wavelet orthonormal basis is 
a family of functions V2^ip{2^x — k)j^ke'^ built by dilations 
and translations of a unique function ijj{x): the analysing 
wavelet. The decomposition thus obtained is this MRA. 
The signal can be reconstructed from this representation 
without any difficulty. We give a minimal background to 
this technique in appendix [Xl 

Now we need to specify the choice of the analysing 
wavelet. As in many problems of image processing, we 
choose the Daubechies wavelet because of its compact sup- 
port. This property is important since it minimizes border 
effects and interactions between patterns of the signal dur- 
ing the filtering process. Moreover, using these wavelets 
also preserves the location of zero-crossing and maxima of 
the signal during the analysis, a property that results in 
mutual suppressive interactions across its different scale 
represe ntations and superior r obustness in noisy environ- 
ments ( Sahiner fc Yaglelll993l) . Thus, the contours of the 
image can be determined efficiently. We understand that 
this property is important in image processing since the 
features of the image are preserved after filtering. For the 
velocity fields we are dealing with, this is also an interest- 
ing point, because we wish to identify flow structures like 
divergences and vortices. 

Finally, wavelet analysis also allows us to determine 
the relevance (or the reality) of flow structures on different 
scales. One may indeed apply the MRA to the velocity 
field and the noise field. Then, for each scale of the flow, 
we can compare the details and the amount of noise to 
see whether the details are relevant or are simply noise 
structures. 



2.4. Curl and divergence fields 

Once an approximation of the velocity field is known, 
it is useful to detect fiow patterns that may be impor- 
tant for the dynamics of the fiuid. As the (measured) 
velocity field is purely two-dimensional, two quantities 
are relevant for enhancing flow structures: the divergence 
D ~ dxVx + dyVy and the z-component of the vorticity 

C = dxVy - dyVx. 

The way derivatives can be computed can be explained 
with a one-dimensional example. Let us consider the ap- 
proximation on scale j of the signal / 



(1) 



where k represents the position of the wavelet. 
Differentiating this expression yields 



dx 



dx 



In the Galerkin method, this formula would be sufficient; 
however, in MRA the derivative of a function approx- 
imated with some resolution has meaning only in the 
same functional space, that is, with the same resolution. 
Therefore, dfj/dx also needs to be projected onto that 
space. Thus we are interested in 

dx 

It therefore turns out that the discrete approximation of 
the derivatives can be easily derived from that of the 
original functions by a simple matrix multiplication. If 
~ {f\^k) discrete approximation of / on scale j 

and S'j!! the discrete approximation of df /dx on the same 
scale, then it follows that 



E 



rk-i 



(2) 



where 



/4>{x — l)-^dx. 
-oo dx 



These numbers can be computed through an algorithm 
described in Beylkin (1992). 

We show the computations of the curl field on a simple 
given velocity field in Fig. [51 namely Vx — ~y, Vy — x. We 
see that, except for the border effects due to the finite size 
of the filter, value 2 is correctly restored. 

2.5. The role of 'holes' 

One of the problems arising when reconstructing the ve- 
locity field comes from the presence of bins without data. 
Such bins produce structures in the divergence and curl 
fields. A simple illustration of the effect is given in Fig. [5] 
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Fig. 5. Isoline plot for the z-component of the curl of the velocity field Vx = —y, Vy = x with some values of v randomly 
set to zero (black dots in the field). To the left, the scafing function 4^ of Daubechies wavelet is used; to the right 
we use 80- The main difference between these two wavelets is the width of their support, twice larger on the right. 
Note the border effects in both figures as well as the patterns introduced by the absence of data and their dependence 
on the support of the wavelet. The large blank areas are at the constant value of 2, as expected; solid lines represent 
isolines of a value different from 2 and dotted lines show negative value isolines. X and Y units are grid points. 



Fig. 6. To the left, the velocity field with a small mesh size (713 km) is complete with 83% of data, while on the right, 
using a larger mesh size of 1223 km, the velocity field is complete to 99%. Patterns of velocities are easily identified 
between the two. We used a time interval of 15 min. 



where the curl of a solid rotation velocity field is plotted. 
This figure shows the importance of the compact support 
of the wavelet in limiting the propagation of errors. 

Let us now consider some real data ta ken from the 
Pic du Midi data set fsee lRoudier et aL 199£l for details). 
Considering the velocity field first, the effects of empty 



bins is not dramatic as may be seen in Fig. [SI velocity 
patterns are indeed not affected much by holes. This is not 
the case for the divergence where we see, in Fig. [71 that 
when the number of empty bins is increased, the patterns 
hardly remain identifiable. Only after a strong filtering 
can one recover similar patterns. 
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Fig. 7. Divergence of the velocity field shown in Fig. [S] As in Fig. [51 on the left we have the more resolved, but less 
complete, sample and on the right the less-resolved, but more complete, sample. In the first row the data have been 
only slightly filtered (they are projected on the space j=l), and one can hardly identify any common structure in the 
two plots. However, we clearly see some common features between the two plots when considering the low frequency 
j=3-component. 



The foregoing example illustrates what may be at the 
origin of the worrying differences found by 'S imon et al.l 
between feature tracking and LCT, namely an in- 
terpolation problem emphasised by considering differen- 
tiated fields. The two methods will converge to similar 
results only when these effects are overcome, which means 
that time or space-averaging is large enough. Indeed, 
empty bins disappear either when the sampling grid is 
coarse enough or when the time resolution is low enough 
(see end of Sect. [2^ . 



3. Discussion and conclusions 

In this paper we have presented the coherent structure 
tracking algorithm aimed at reconstructing the horizontal 



velocity field on the solar surface from the granule mo- 
tions. 

We first discussed the role of segmentation and de- 
scribed a way to generalise Strous algorithm. We have 
shown that a segmentation that follows the shape of gran- 
ules more closely is more sensitive to the splitting of large 
granules and that, as far as velocity measurements are 
concerned, the Strous criterion remains the most efficient. 
We then showed that the reconstruction of the velocity 
field is a delicate process because it is subject to many 
constraints. Indeed, granules do not sample the field of 
view uniformly, and reconstruction of the velocity field, 
along with its derivatives like divergence or curl, requires 
some interpolation. There are many ways of performing 
such an operation; however, classical methods, like poly- 
nomial interpolation, would propagate errors and noise 
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everywhere. We thus selected a method based on MRA 
which projects data onto Daubechies wavelets. The finite 
support of these functions hmits the effects of noise and 
error propagation: sides and regions lacking in data have 
a limited influence. Moreover, the signal is decomposed at 
different scales through a filtering process. At each step 
the filtered field and the remaining details can be viewed 
and compared. This representation is particularly relevant 
for turbulent flows, since the relation between amplitude 
and scale is of crucial importance for constraining models 
of these flows. 

We also related the minimum size of the velocity mesh 
grid to the time resolution. We thus found that, typically, 
one granule trajectory occupies a "volume" of 1200Mm^s. 
When the "space-time" resolution does not reach this 
limit, many granules contribute to the velocity in one mesh 
point. Their mean velocity is considered as the true lo- 
cal velocity, but local fluctuations around this mean gives 
some information on the local strength of convection. 

We did not discuss here the influence of the Earth's 
atmospheric turbulence on the determination of the ve- 
locity fields and derivatives. Surely this is an important 
point: typical atmospheric distortion of images induces, 
in good seeing conditions, feature motions of 0'.'13. For 
a long-lived structure, say 10 min, this means an uncer- 
tainty of the velocity as high as 220 m/s. This is quite 
large compared to a typical velocity of 600 m/s. Hence, 
atmospheric noise is a non negligible part of the data and 
a careful determination of its influence is needed. This is 
the subject of the companion pa per to which we refer the 
reader (see iTkaczuk et al~ 2007 ). 

Finally, although this has not been tested yet, we think 
that the CST algorithm can be fruitfully used to track the 
magnetic features of the photosphere, like network bright 
points or even determine velocities of features in the solar 
atmosphere. 



Appendix A: Fundamentals of multi-resolution 
analysis 

We give in this appendix the basic background of MRA 
and refer the r eader to textbooks for a more co mplete pre- 
sentation fe.g. lDaubechieslll992l: iMallatI [l999l ) . An MRA 
is a sequence {Vj}-^-^ of closed subspaces of square- 
integrable functions, L^(IR), such that the five following 
properties are satisfied, Vj G 2: 



1. V,CV,+^. 

Vj can be interpreted as the set of all possible signal 
approximations at the resolution 2^ (a resolution r is 
defined by the size 1/r of the smallest detail). Thus, 
if the smallest detail in Vq has size 1, it is possible to 
read on Vj details of size 2^^ . 

It follows from this property that the approximation of 
the signal at resolution 2^~^^ contains all the necessary 
information for determining the same signal at a lower 
resolution 2^ (plus additional details). 



2. (J Vj is dense in L^(JR) and Q Vj = {0}. 

In other words, when j increases, the approximated 
signal converges to the original signal. Conversely, if 
j (the resolution) decreases, the approximated signal 
converges to zero (it contains less and less informa- 
tion). 

3. fix) e Vj^fi2x) e Vj+i. 

This property defines 2 as the rate of scaling change 
(the ratio of two successive resolution values) . 



4. fix) e Vj^fix-2-^k) e Vj 



Vfc e z. 



This property characterizes the invariance under dis- 
crete translations: when the signal is translated by 
a length proportional to 2~^, the approximations are 
translated by the same amount and no information is 
lost in the translation. 
5. A function exists in Vq such that {(/)(a; — fc)}fcg2 
is an orthonormal basis of Vq. Hence, the family 
{2^/20(2^2; - fc)}feg2 is an orthonormal basis of Vj. 
Function (jj is called the scaling function of the multi- 
resolution representation. 



Now, let us define Wj as the orthogonal complement of Vj 
in Vj+i (it contains the additional details that are in V^+i 
and not in Vj). There exists a function (the wavelet) 
such that {V'(2^^fc)}feez is an orthonormal basis of Wq and 
{2^^^ip{2^ X — k)}j ke:E is an orthonormal basis of ^^(IR). 

One studies a signal / of L^(]R) by projecting it or- 
thogonally on the collection of Vj and Wj . This procedure 
can be carried out according to the pyramidal algorithm 
presented below. 

First, let us define the two filters h and g that can 
be deduced from the MRA. This analysis allows us to de- 
termine a function h, which is the impulse response of 
some 27r-periodic low-pass filter H defined with the scal- 
ing function: Hiuj) — (j>i2uj)/(j>ioj). On the other hand, 
one defines function g by Giuj) — il>i2uj) / (j)iuj), g being 
the impulse response of the 27r-periodic high-pass filter 
G. The filters H and G are quadratic mirror filters and 
are finked by the relation Giu) = e~''^if(a; -I- tt), giving 
gin) = (— l)^~"/i(l — n) for the impulse responses. Then, 
the pyramidal architecture for computing the wavelet rep- 
resentation can be easily written as: 



suppose that /(x^) belongs to the Vq ifixi) = foixi) 
approximation of the signal at resolution 1) and de- 
compose fixi) onto V-i and W-i; 
the decomposition onto V-i consists in a convolution 
by the filter ft,, such that /i(n) = /i(— n), and a decima- 
tion (i.e. only one out of every two sample is retained); 
we obtain the 7V/2— sampled so-called approximation 
at resolution 1/2 equal to 



f-i{n) 



N 

E 



ft(2fc- n)/o(n) e VI 
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— the decomposition onto W^i consists in a convolution 
by filter g (such that g{n) — g(—n)) and a decimation; 
we obtain in the same way 



N 



this is the detail at resolution 1/2, that is to say the 
"difference in information" between fo{n) and 
it also has N/2 samples. 

By repeating the same sequence, we obtain the approxi- 
mation and the detail at resolution 1/2^ : 

N/2 

f-2{n) = YJ^{2k^n)f^^{n) € y_2 

n=l 

and 

N/2 

d_2(n) = ^3(2A;-n)/_i(n) £ W.2, 

n=l 

and so on. 

After a number J of iterations to be defined by the 
problem, we have decomposed /o into d-i, d_2, • ■ • , dj, 
and fj. 

Let us finally mention that reconstruction of fo{n) 
from the details and the last approximation is just as easy 
and has appreciable quality. One has to iterate (starting 
from J = J) : 



M 



M 



f,{n) = 2(J2 Hn - 2fc)/,_i(fc) + ^ g(n - 2fc)d,.i(fc) 



fc=i 



fc=i 



fj-i and dj-i being sampled in M points. 

The MRA can be generalised to two dimensions for 
image- processing applications. We can define a sequence 
of multi-resolution vector sp aces and the approximations 
of a signal /(x, y) € L'^(5R^) (|Mallatlll989[) . Since the im- 
age under study is bounded, we choose (j) with a compact 
support (i.e. vanishes outside a finite region). In that case, 
filters h and g have only fi nitely many coeffic ients that sat- 
isfy the MRA conditions ( Daubechied 199^ ). Functions (j) 
and -0 become more regular as the number of coefficients 
n increases (the case of n = 1 correspo nds to the discon - 
tinuous Haar basis). It has been proved (|Daubechieslll992h 
that the regularity of (j) and increases linearly with n. 
By choosing n = 8, we obtain a good compromise between 
regularity and support width. 



Mallat, S. 1999, A guided tour of signal processing 

(Academic Press) 
Meyer, Y. 1987, in wavelets, time-frequency methods and 

phase space, ed. J. M. Combes, A. Grossmann, & 

P. Tchamitchian (Springer Verlag), 21-37 
MuUer, R., Auffret, H., Roudier, T., et al. 1992, Nature, 

356, 322 

November, L. J. & Simon, G. W. 1988, ApJ, 333, 427 
Potts, H. E., Barrett, R. K., & Diver, D. A. 2003, 

Sol. Phys., 217, 69 
Rieutord, M., Roudier, T., Ludwig, H.-G., Nordlund, A., 

& Stein, R. 2001, A & A, 377, L14 
Roudier, T., Rieutord, M., Malherbe, J., & Vigneau, J. 

1999, A & A, 349, 301 
Sahiner, B. & Yagle, A. E. 1993, IEEE Trans, on Signal 

Proc, 41, 3579 
Simon, G. W., Brandt, P. N., November, L., Shine, R., 

& Strous, L. 1995, in Proc. of 4th SOHO Workshop: 

Helioseismology (ESA SP-376), 223 
Soille, P. 1999, Morphological Image Analysis (Berlin: 

Springer) 

Strous, L. 1995, in Proc. of 4th SOHO Workshop: 

Helioseismology (ESA SP-376), 213 
Tkaczuk, R., Rieutord, M., Meunier, N., & Roudier, T. 

2007, A & A, to appear, 1 
Vincent, L. & Soille, P. 1991, IEEE Transactions on 

Pattern Analysis and Machine intelligence, 13, 583 



References 

Adrian, R. J. 2005, Experiments in Fluids, 39, 159 
Beylkin, G. 1992, SIAM J. Numer. Anal., 6, 1716 
Bovelet, B. & Wiehr, E. 2001, Solar Phys., 201, 13 
Daubechies, I. 1992, Ten lectures on wavelets (SIAM) 
Mallat, S. 1989, IEEE Trans, on Pat. An. and Mach. Int., 
11, 674 



